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1.  Introduction 


Future  military  engagements  will  require  weapons  systems  exhibiting  improved  range  and 
accuracy.  One  of  the  technologies  under  investigation  to  achieve  these  goals  is  the 
electrothermal-chemical  (ETC)  propulsion  concept  shown  schematically  in  figure  1.  In  the  ETC 
gun,  energy  that  is  stored  either  in  batteries  or  in  a  rotating  device  is  converted  on  demand  into 
an  electrically  generated  plasma  (resulting  from  the  ablation  of  polyethylene  material  in  a 
capillary)  that  is  injected  into  the  chamber  in  a  howitzer  or  gun.  This  plasma  energy  is  used  to 
ignite  the  chemical  propulsion  charge  (i.e.,  solid  propellant)  as  well  as  enhance  gun  performance 
by  taking  advantage  of  a  number  of  unique  plasma  characteristics.  For  example,  a  low-density 
plasma  jet  can  efficiently  ignite  charges  of  high-loading  density,  can  control  propellant  mass 
generation  rates  (1),  can  reduce  propellant  charge  temperature  sensitivity  (i.e.,  the  variation  of 
gun  performance  with  changing  ambient  temperature  (2,  3)  and  can  shorten  ignition  delay 
(i.e.,  the  time  interval  between  firing  of  the  igniter  and  ignition  of  the  propellant  (4).  Plasma 
igniters  eliminate  the  conventional  chemical  igniter  and  thus  can  enhance  the  safety  aspects  of 
the  overall  gun  propulsion  system. 


Figure  1.  ETC  gun  concept. 


Research  has  been  carried  out  on  many  aspects  of  plasmas  related  to  uses  such  as  solid  propellant 
igniters  (5,  6).  Because  the  plasma  is  at  a  temperature  (typically  >10,000  K)  that  is  considerably 
higher  than  conventional  chemical  igniters  (3000  K),  the  radiation  properties  of  the  plasma  have 
also  been  considered.  The  high  plasma  temperature  leads  to  radiation  effects  nearly  100  times 
greater  than  that  of  chemical  igniters  (i.e.,  a  T4  effect)  (5).  Such  radiation  could  lead  to 
significantly  different  temperature  profiles  within  the  propellant,  causing  changes  in  combustion 
rates.  Plasma  has  a  much  lower  density  than  the  gases  generated  by  a  chemical  igniter.  This  will 
impact  the  convective  heat  transfer  to  the  propellant,  as  the  plasma  moves  through  the  grains  as 
well  as  the  velocity  and  mode  of  flamespreading  within  a  propellant  bed.  It  has  been  suggested 
that  energy  transport  by  convection  may  be  as  important  as  radiation  transport  in 
plasma-propellant  interactions  (PPI)  (7,  8). 
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All  of  the  effects  described  previously  can  lead  to  significant  changes  in  ballistic  behavior  and 
useful  improvements  in  gun  performance,  but  an  understanding  of  the  underlying  physical 
mechanisms  is  necessary  to  achieve  this  goal.  To  this  end,  the  U.S.  Army  Research  Laboratory 
(ARL)  has  begun  a  comprehensive  study  on  the  interaction  of  the  plasma  efflux  from  an  ETC 
igniter  with  solid  propellant  grains.  The  goal  of  this  work  is  to  elucidate  the  relevant  physical, 
mechanical,  and  chemical  mechanisms  that  underlie  the  observed  ballistic  effects.  Various 
aspects  of  the  experimental  program  for  PPI  are  described  elsewhere  (9, 12).  The  first  phase  of 
the  modeling  effort  in  support  of  this  PPI  project  involves  a  time-accurate  computational  fluid 
dynamics  (CFD)  code  that  has  been  written  to  include  high-temperature  thermodynamics, 
variable  specific  heats  and  transport  properties  (viscosity  and  thermal  conductivity),  and 
finite-rate  (nonequilibrium)  chemical  kinetics  (the  mechanism  is  described  in  reference  (10).  A 
separate  capillary  model  described  by  McQuaid  and  Nusca  (11)  supplies  boundary  conditions  for 
the  CFD  code  in  terms  of  the  physical  and  chemical  properties  of  the  plasma  capillary  efflux. 
Validation  of  the  capillary  model  and  the  CFD  code,  including  coupled  chemistry,  is  conducted 
by  simulating  a  series  of  experiments  wherein  a  plasma  jet  is  generated  from  a  plasma  cartridge; 
pressures  in  the  resulting  unsteady  flowfield  are  measured  using  probes  mounted  on  a  plate  held 
normal  to  the  efflux. 


2.  Experimental  Efforts  for  PPI 


White  et  al.  (8)  and  Williams  and  White  (12)  documented  the  plasma-propellant  interaction 
experiments  conducted  at  ARL.  The  general  view  of  the  plasma  generator  is  shown  in 
figure  2.  The  plasma  capillary  consists  of  a  standard  polyethylene  liner  with  dimensions 
102  mm  in  length  and  6.35  mm  in  diameter.  The  nozzle  consists  of  a  circular  tube  of  2.9  cm  in 
length  and  1.3  cm  in  diameter  that  is  directly  adjacent  to  the  end  of  the  plasma  capillary  followed 
by  another  circular  tube  of  1.9  cm  in  length  and  3.2  cm  in  diameter.  As  a  result,  the  plasma 
efflux  generated  by  the  capillary  is  expanded  through  a  step-configuration  nozzle  before 
emerging  into  the  open  air.  A  pressure  probe  is  positioned  in  the  path  of  the  plasma  efflux 
(figure  2).  To  protect  this  probe,  the  plasma  capillary  is  mounted  behind  a  steel  plate.  The 
pressure  probe  was  placed  15  cm  away  from  the  surface  of  the  plate  and  aligned  with  the 
centerline  of  the  cartridge.  Care  was  exercised  in  shielding  the  gage  from  direct  interaction  with 
the  plasma  as  were  the  cables  and  amplifiers  (strong  evidence  of  charged  particle  impact  on  the 
unshielded  gage  was  found  [5]).  Even  so,  post-firing  examination  of  the  pressure  gage  indicated 
significant  amounts  of  particulate  in  the  plasma  effluent  (the  plasma  is  generated  by  ablation) 
that  contributed  to  noticeable  variability  in  pressure  histories  from  shot  to  shot.  Currently, 

R.  A.  Beyer  of  ARL  is  conducting  a  new  set  of  experiments  in  which  the  plasma  efflux  is 
injected  into  a  closed  chamber  (13). 
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Figure  2.  Schematic  of  ARL  open-air  experiment  (also  showing  computational  domain). 


Litzinger  et  al.  (14)  of  Penn  State  University  (PSU)  have  designed  and  operated  two  PPI 
experiments.  The  open-air  experiment  is  shown  schematically  in  figure  3.  The  polyethylene 
capillary  is  typically  26  mm  in  length  and  3  mm  in  diameter  (d).  An  extension  tube  26  mm  in 
length  is  placed  at  the  capillary  exit,  which  guides  the  plasma  efflux  into  the  open-air.  The 
capillary  and  the  extension  tube  are  mounted  within  a  solid  housing  (not  shown  in  figure  3) 
making  the  setup  similar  to  that  shown  in  figure  2  (i.e.,  the  protective  plate  mounted  near  the 
capillary  exit).  An  instrumented  plate  is  placed  at  some  distance  (L)  from  the  plasma  device, 
typically  19  mm.  Pressure  probes  are  mounted  on  the  plate  with  a  spacing  of  9.53  mm  (3/8  in). 


Anode 

(Tungsten  alloy) 


d=l/8"  or  1/4" :  bore  diameter 
a=3/8"(9. 525mm):  seperationofPTs 


Figure  3.  Schematic  of  the  PSU  experimental  setup. 
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3.  Previous  Modeling  Efforts 


The  initial  modeling  effort  for  open-air  plasma  discharges  was  conducted  at  ARL  from  1997  to 
1998  (15-16).  This  computational  study  was  based  on  the  National  Research  Laboratory  (NRL) 
CFD  code  FAST3D  (17)  and  sponsored  by  the  Department  of  Defense  (DOD)  Common 
High-Performance  Computing  (HPC)  Software  Support  Initiative  (CHS  SI)  program.  The  ARL 
open-air  plasma  experiment,  shown  schematically  in  figure  2,  was  chosen  as  a  test  case  for  the 
FAST3D  CFD  code  because  the  plasma  jet  is  of  very  high  velocity  and  temperature,  but  very 
low  density.  The  plasma  jet  was  assumed  to  be  an  ideal,  neutral,  single-component  gas  with 
density,  temperature,  ratio  of  specific  heats  (y),  and  velocity  characteristics  identical  to  the 
plasma  discharge.  The  multicomponent  nature  of  the  plasma,  consisting  of  various  neutral 
species,  ions,  and  electrons,  and  its  chemical  reactivity  were  not  considered  in  those  studies.  The 
augmentation  of  plasma-air  mixing,  caused  by  viscosity  and  turbulence,  was  also  excluded  in 
these  studies.  The  time-dependent  gas  properties,  along  with  jet  velocity  and  y  were  determined 
using  a  plasma  capillary  code  developed  at  ARL  by  J.  Powell  (18).  Overall,  these  CFD 
simulations  captured  the  major  features  of  the  plasma  jet  as  photographed  by  White  et  al.  (5). 
Alternating  bright-dark  structures  in  the  photos  were  identified  as  the  shock  structures  in  the  jet. 
A  precursor  shock  caused  by  the  high-velocity/low-density  plasma  effluent  (a  feature  not  easily 
photographed)  was  captured.  Comparisons  between  computed  and  measured  pressures  in  the  jet 
were  encouraging  but  pointed  to  the  need  for  a  variable  y,  multicomponent,  reacting  flow  model. 

More  advanced  modeling  was  conducted  by  Nusca  et  al.  (19)  in  1999-2000  and  is  continued  in 
the  present  effort  (see  section  4  for  a  description  of  the  model).  Figures  4  and  5  show  some 
typical  results  from  that  effort.  The  ARL  open-air  plasma  jet  experiment  was  used  for  code 
validation.  Figure  4  demonstrates  the  model’s  utility  in  matching  the  time-of-arrival  of  the 
precursor  shock  (shown  just  passing  over  the  probe  tip  in  figure  5)  and  the  stagnation  pressure  of 
~7  MPa  (pressure  values  in  figure  5  are  nondimensional  so  that  a  value  of  40  corresponds  to  a 
pressure  of  ~4  MPa).  Although  these  comparisons  are  good,  it  was  felt  that  a  better  agreement 
between  predicted  and  measured  pressure  values  was  required.  This  report  discusses  a 
comparison  between  this  model  and  the  Penn  State  open-air/plate  experiment 
(figure  3). 


4.  Multispecies  Reacting  Flow  CFD  Code 


The  high-temperature,  nonideal  chemically  reacting  gas  flowfield  within  the  capillary  efflux  jet 
is  numerically  simulated  using  CFD.  The  NSRG2  code  (20-22)  solves  the  2-D/axisymmetric, 
unsteady  real-gas  Navier-Stokes  equations  including  submodels  that  represent  finite-rate 
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Figure  4.  Measured/computed  pressure  at  the  probe  tip. 


Figure  5.  Computed  pressure  contours  at  probe. 

(nonequilibrium)  chemical  reactions,  multispecies  diffusion,  as  well  as  variable  specific  heats, 
viscosity,  and  thermal  conductivity.  These  partial  differential  equations  are  cast  in  conservation 
form  and  converted  to  algebraic  equations  using  a  finite-volume  formulation.  Solution  takes 
place  on  a  mesh  of  nodes  distributed  in  a  zonal  fashion  throughout  the  flowfield  such  that  sharp 
geometric  comers  and  other  details  are  accurately  represented.  The  conservation  law  form  of  the 
equations  ensures  that  the  end  states  of  regions  of  discontinuity  (e.g.,  shocks)  are  physically 
correct  even  when  smeared  over  a  few  computational  cells.  The  Navier- Stokes  equations  for 
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2-D/axisymmetric  (x,  y  coordinates),  reacting  (N  species)  and  unsteady  (time,  t)  flow  are  written 
in  nondimensional  form  (20).  The  dependent  variables  are  density  (p),  velocity  ( V  and 
components  u,  v),  energy  (e),  and  species  mass  fraction  (cl).  Pressure  and  temperature  are  p  and 
T,  whereas  L  is  a  characteristic  lengthscale  for  the  flowfield. 
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where  ayy  and  qy  are  similarly  written  and  E  =  V  lh,  =  hr  A,  A  =  pLe/RePr,  Le  =  Pr/Sc, 

Pr  =  pCp/fc,  Sc  =  pi pD,  and  Re  =  pVL/p.  The  Navier-Stokes  equations  (equation  1)  are  written  in 
integral  form  and  then  re-expressed  in  a  semi-discrete  fashion  using  a  finite-volume 
discretization  technique.  The  numerical  computations  proceed  by  solving  the  semi-discrete 
equation  on  each  computational  cell  using  central  and  upwind  numerical  differences  along  with 
flux  limiting.  Once  properly  discretized,  the  resulting  set  of  algebraic  equations  is  solved  in  a 
coupled  manner  in  time  using  an  explicit  time-accurate  method.  The  numerical  time-step  is 
computed  using  the  CFL  condition.  A  separate  chemical  time-step  is  computed  as  well.  The 
final  time-step  is  based  on  the  smallest  of  these.  Nusca  (20-22)  provides  a  more  complete 
description  of  the  numerical  scheme. 
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The  diffusion  coefficient,  D,  is  understood  to  be  the  multicomponent  coefficient,  Dim  (specie  i 
into  the  mixture),  which  is  related  to  the  binary  diffusion  coefficient,  Ay  (specie  i  diffusing  into 
specie/), 


A .  = 


1  -X 


y  x  ./d 

Lu j  j  v 


Dn  = 


0.00266  T/ 2 


pM< 


(2) 


where  X,  is  the  mole  fraction  and  M,  is  the  molecular  weight  of  specie  i,  Q d  is  the  collision 
integral  and  o(/  is  the  collision  diameter  (available  in  the  literature).  The  thermal  diffusion  ratio, 
kT ,  is  assumed  to  be  1 .0.  Because  no  provisions  have  been  made  in  the  conservation  equations 
(equation  1)  for  flows  with  electric  currents  at  this  time,  the  flowfield  was  rendered  electrically 
neutral  by  setting  the  diffusion  velocity  of  the  electrons  equal  to  the  average  diffusion  velocity  of 
the  ions — the  diffusion  coefficient  for  the  electrons  is  then  computed  from  that  of  the  ions. 

Future  works  will  be  aimed  at  relaxing  this  assumption  by  allowing  for  charge  separation  effects 
and  including  appropriate  source  terms  in  the  momentum  and  energy  equations. 


We  realize  that  the  plasma  gas  does  not  necessarily  behave  as  a  perfect  fluid.  Indeed,  the  ionized 
gas  is  usually  characterized  as  rarefied  and  one  in  which  Coulomb  interactions  between  charged 
particles  create  significant  departures  from  the  perfect  gas  behavior.  However,  for  “weakly” 
imperfect  gases,  one  can  prescribe  terms  that  account  for  Coulomb  interactions  as  corrections  to 
the  classic  equations.  Such  corrections  are  described  in  reference  (23).  These  corrections  can  be 
expressed  in  the  pressure-density-temperature  relation  using 


P  =  PRTt,17 

1  =  1  M, 


P 


correct 


(3) 


where  z,  is  the  charge  of  the  i-th  species  and  the  constant  K  =  2.4703  x  1010  for  p  in  N/m2. 
Whereas  the  energy,  density,  and  mixture  properties  (mass  fractions  and  molecular  weight)  are 
determined  from  the  conservation  equations,  the  temperature,  T,  to  be  used  for  equation  3  is 
decoded  from  the  energy  using  the  local  density  and  mixture  properties  and  by  iteratively  solving 


e 


h(T,Ci ) 


RT 


^correct  (P  ’  ^  ’  Ci) 


(4) 


where  the  enthalpy,  h,  is  determined  using  the  NASA-Lewis  database  (24)  and  the  constant 
K  =  2.4703  x  107  for  e  in  joules  per  kilogram.  It  has  since  been  shown  (23)  that  the  property 
determinations  obtained  from  the  National  Aeronautics  and  Space  Administration 
(NASA)-Lewis  database  (24)  without  modification,  such  as  the  enthalpy  and  specific  heats,  are 
reasonably  accurate  and  that  the  correction  terms  need  not  be  employed.  Thus,  use  of  such 
corrections  will  not  be  continued  in  this  effort.  Rather,  the  mixture  temperature  and  pressure  are 
computed  using 
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2  2 

U  +  V 


N 


N 


dc. 


2>,  W°fi  \;c,=  X  V,+2>,  [  ;  p  =  p Z  *  r£^-,  (5) 


!  =  1 


=  1 


=  1 


where  Z  is  the  compressibility  factor  for  the  mixture  (Z  <  1).  The  use  of  one  gas  temperature  for 
all  species  implies  thermal  equilibrium  even  though  chemical  nonequilibrium  is  utilized  in  the 
CFD  model. 


The  species  viscosity  and  thermal  conductivity  coefficients  are  given  by  polynomials  with 
coefficients  A,  B,  C,  and  D  given  in  reference  (24)  (coefficients  are  different  for  each  variable 
and  each  species).  The  mixture  values  of  viscosity  and  thermal  conductivity  are  computed  using 
the  values  for  each  species  and  empirical  mixture  rules  discussed  in  reference  (20).  For  each 
species,  i, 


Hb)  =  \.i  HT)  +  +  jf 


+  D„ 


ln(x)  -  AKi  ln(T)  +  +DKi  .  (6) 


A  chemical  source  term  appears  in  equation  1  for  each  species,  i,  denoted,  vv,.  Chemical  reactions 
can  be  expressed  in  a  general  reaction  rate  equation  with  stoichiometric  coefficients  for  species 
in  these  reactions.  The  chemical  source  terms  that  appear  in  the  H  array  for  equation  1  (i.e.,  >v) 
are  computed  using  this  general  reaction  rate  equation: 


dC  _  wj 
dt  M  t 


(A  -  v-;j  k  n c 


N 


k  n  cr- 


(7) 


where  the  C  above  represents  concentration  that  is  raised  to  powers  of  either  the  (v-prime)  or 
product  coefficient  (v-double-prime).  The  forward  reaction  rate  is 


kf  =  CT"  exp  (-E/kT) ,  (8) 

where  the  rate  data  C,  E,  and  n  are  given  in  table  1 .  The  backward  rate  is  computed  using  the 
equilibrium  constant  for  each  reaction 

K  =kf!Kp  ;  Kp  =  exp  (-A%r),  (9) 

where  Kp  is  the  equilibrium  constant  for  a  particular  reaction  computed  from  the  change  in  Gibbs 
energy  for  that  reaction  (see  reference  [20]  for  details).  Gibbs  energy  for  each  species  is 
computed  using  the  NASA-Lewis  tables  (see  reference  [24]). 

Where  equilibrium  chemistry  is  assumed  inside  the  capillary,  resulting  in  a  mixture  primarily 
composed  of  electrons,  C,  C+,  H,  and  H+  emerging  from  the  orifice,  a  chemical  kinetics 
mechanism  was  developed  to  describe  the  finite -rate  chemistry  that  takes  place  outside  the 
orifice  of  the  capillary  as  the  emerging  jet  mixes  with  air  (10).  An  important  simplifying 
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assumption  made  was  that  it  is  sufficient  to  consider  only  mixtures  which  had  a  concentration  of 
O2  much  greater  than  that  of  the  plasma  constituents;  that  is,  the  mixture  considered  was  fuel 
lean  in  the  combustion  sense,  the  C  and  H  containing  species  mentioned  previously  being  fuels. 
We  hope  to  subject  this  assumption  to  later  scrutiny  and,  possibly,  revision  as  more  complex 
models  are  developed  in  the  future.  For  initial  modeling,  this  assumption  was  made  because 
fuel-rich  chemistry,  as  might  occur  in  the  plasma  jet,  is  known  to  be  extremely  complex.  The 
notion  that  the  important  mixtures  are  lean  may  be  reasonable  because  typical  masses  of  plasma 
formed  during  an  ETC  gun  cycle  are  much  smaller  than  that  of  the  air  in  the  chamber  ullage  (25). 

The  chemical  mechanism  was  developed  in  two  stages.  First,  a  detailed  mechanism  consisting  of 
224  reactions  involving  5 1  species  was  compiled.  The  core  of  this  mechanism  was  derived  from 
a  chemical  mechanism  developed  over  the  last  10  years  for  modeling  the  dark  zones  of 
combusting  solid  propellants  (26).  To  this  were  added  reactions  of  our  choosing  that  are 
believed  to  be  the  most  important  involving  the  charged  species  and  C  atoms  present  in  the 
plasma  (these  species  are  not  important  in  propellant  dark  zones).  The  addition  of  reactions 
beyond  those  included  in  the  dark  zone  mechanism  were  largely  guided  by  thermodynamic 
considerations  and  literature  review.  Most  of  the  rate  constants  could  be  obtained  from  the 
literature.  However,  some  of  the  reactions  thought  to  be  important  have  never,  to  our 
knowledge,  been  studied.  Rate  constants  for  these  had  to  be  estimated.  Fortunately,  the 
estimated  reactions  are  all  highly  exothermic.  As  a  result,  it  is  felt  the  estimates  are  likely  to  be 
quite  good  (within  a  factor  of  five).  Also,  because  the  reactions  are  expected  to  be  extremely 
fast,  it  is  unlikely  modeling  results  for  gross  properties  will  be  very  sensitive  to  the  rate 
coefficients  assumed. 

The  second  stage  in  chemical-mechanism  development  was  reduction  of  the  detailed  chemical 
mechanism  to  a  smaller  skeletal  mechanism.  The  skeletal  mechanism  retains  all  the  elementary 
reactions  found  to  be  important  in  describing  the  chemical  evolution  of  a  lean  plasma/air  mixture 
across  a  wide  range  of  flowfield  temperatures  and  pressures.  This  smaller  mechanism  consists  of 
57  reactions.  Table  1  lists  the  reactions  and  rate  data  (note  that  the  units  for  C  are  mole,  second, 
centimeter  cubed,  and  Kelvin,  EJk  is  in  units  of  K;  n  is  nondimensional,  whereas  k  is  the 
Boltzmann  constant). 

The  rate  for  reaction  5  is  of  the  additive  type  wherein  two  rates  are  computed  (for  reactions  5  and 
5a)  and  added  together  for  the  final  rate.  The  rates  for  reactions  30,31,  and  55  are  of  the 
“falloff  ’  type,  wherein  the  final  rate  is  the  product  of  three  factors:  two  separate  rates  (e.g., 
reactions  3 1  and  31a)  and  a  function  based  on  the  local  flowfield  temperature  and  mixture,  F 
(T,  mixture).  See  reference  (10)  for  further  details.  The  rate  data  for  reactions  5a,  30a,  31a,  and 
55a  are  given  in  table  2.  Certain  reactions  involve  a  “third-body”  or  a  “collision  partner,”  which 
is  denoted  M.  The  species  M  can  stand  for  any  of  the  species  being  considered  in  the 
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Chemical  kinetics  mechanism  (10). 

"T  Reaction  T 


Comments 


1 

H+  +  O2  —  H  +  O2+ 

1.80E+  15 

0 

0 

— 

2 

C+  +  02=C0+  +  0 

3.00E+  14 

0 

0 

— 

3 

C+  +  02  =  0+  +  CO 

3.00E+  14 

0 

0 

— 

4 

H+  +  C  =  H  +  C+ 

1.80E+  15 

0 

0 

— 

5 

0+  +  O2  —  0  +  O2+ 

8.27E+  14 

-0.77 

0 

See  additive 
rates;  5a 

6 

0+  +  N2  =  NO+  +  N 

6.00E+  11 

0 

0 

— 

7 

CO+  +  0  =  0+  +  CO 

8.40E  +  13 

0 

0 

— 

8 

e  +  O2+  =  20 

4.90E+  18 

-0.65 

0 

— 

9 

e  +  CO+  =  C  +  0 

1.86E+  18 

-0.48 

0 

— 

10 

e  +  NO+  =  N  +  0 

3.10E+  19 

-0.85 

0 

— 

11 

C  +  02  =  CO  +  0 

5.80E+13 

0 

289.88 

— 

12 

02  +  H  =  0  +  OH 

3.52E+  16 

-0.7 

8590.84 

— 

13 

N  +  NO  =  N2  +  0 

3.27E+  12 

0.3 

0 

— 

14 

0  +  NO  =  N  +  02 

3.80E  +  09 

1 

20822.85 

— 

15 

NO  +  H  =  N  +  OH 

1.70E+  14 

0 

24559.64 

— 

16 

CO  +  OH  =  C02  +  H 

1.51E  +  07 

1.3 

-381.48 

— 

17 

H  +  H02  =  20H 

1.69E+  14 

0 

439.86 

— 

18 

H  +  H02  =  H2  +  02 

6.63E+  13 

0 

1069.95 

— 

19 

H  +  02  +  M  =  H02  +  M 

3.61E+  17 

-0.72 

0 

— 

20 

20H  =  0  +  H20 

6.00E  +  08 

1.3 

0 

— 

21 

OH  +  H2  =  H20  +  H 

2.16E  +  08 

1.5 

1726.22 

— 

22 

H  +  OH  +  M  =  H20  +  M 

1.60E  +  22 

-2 

0 

— 

23 

2H  +  M  =  H2  +  M 

1.00E+  18 

-1 

0 

— 

24 

2H  +  H2  =  2H2 

9.20E  +  16 

-0.6 

0 

— 

25 

2H  +  H20  =  H2  +  H20 

6.00E  +  19 

-1.25 

0 

— 

26 

2H  +  C02  =  H2  +  C02 

5.49E  +  20 

-2 

0 

— 

27 

OH  +  H02  =  H20  +  02 

7.50E+  12 

0 

0 

— 

28 

NO  +  H02  =  N02  +  OH 

2.11E+  12 

0 

-241.07 

— 

29 

no2  +  o  =  no  +  o2 

3.90E+  12 

0 

-119.78 

— 

30 

N02  +  M  =  NO  +  0 

7.60E+  18 

-1.27 

36884.75 

See  falloff 
rates;  30a 

31 

NO  +  OH  +  M  =  HN02  +  M 

1.99E+  12 

-0.05 

-362.86 

See  falloff 
rates;  31a 

32 

HN02  +  OH  =  H20  +  N02 

1.27E+  10 

1 

67.94 

— 

33 

co  +  no2  =  no  +  co2 

9.04E+  13 

0 

17000.5 

— 

34 

N02  +  H  =  NO  +  OH 

1.30E+  14 

0 

181.68 

— 

35 

H02  +  CO  =  C02  +  OH 

5.80E+  13 

0 

11542.02 

— 

36 

H202  +  M  =  20H  +  M 

1.30E+  17 

0 

22898.84 

— 

37 

2HO2  =  H2O2  +  O2 

1.80E+  12 

0 

0 

— 

38 

H202  +  OH  =  H20  +  H02 

1.75E+  12 

0 

160.04 

— 

39 

0  +  H02  =  02  +  OH 

1.40E+  13 

0 

540.01 

— 

40 

20  +  M  =  02  +  M 

1.89E+  13 

0 

-899.85 

— 

Table  1.  Chemical  kinetics  mechanism  (10)  (continued). 


No. 

Reaction 

c 

n 

E/k 

Comments 

41 

NO  +  M  =  N  +  0  +  M 

1.40E+  15 

0 

74700.55 

— 

42 

N2  +  M  =  2N  +  M 

3.71E  +  21 

-1.6 

113236.0 

3 

— 

43 

NH  +  M  =  N  +  H  +  M 

2.65E+  14 

0 

38002.01 

— 

44 

NH  +  0  =  NO  +  H 

5.50E+  13 

0 

0 

— 

45 

NH  +  0  =  N  +  OH 

3.72E+  13 

0 

0 

— 

46 

C  +  OH  =  CO  +  H 

5.00E+  13 

0 

0 

— 

47 

0  +  H2  =  OH  +  H 

5.06E  +  04 

2.67 

3165.58 

— 

48 

CN  +  N  =  C  +  N2 

1.04E+  15 

-0.5 

0 

— 

49 

CN  +  0  =  CO  +  N 

2.05E+  13 

0 

209.86 

— 

50 

H  +  0  +  M  =  OH  +  M 

6.20E+  16 

-0.6 

0 

— 

51 

H02  +  03  =  OH  +  202 

8.40E  +  09 

0 

598.89 

— 

52 

0  +  O3  =  202 

4.82E+  12 

0 

2058.38 

— 

53 

H  +  03  =  OH  +  02 

8.43E+  13 

0 

478.11 

— 

54 

NO  +  O3  —  N02  +  02 

1.08E+  12 

0 

1368.9 

— 

55 

0  +  02  +  M  =  O3  +  M 

3.49E+  12 

0 

231.5 

See  falloff 
rates;  55a 

56 

0+  +  NO  =  NO+  +  0 

1.00E+  12 

0 

0 

— 

57 

02+  +  NO  =  02  +  NO+ 

3.00E+  14 

0 

0 

— 

Table  2.  Rate  data  for  special  reactions. 


No. 

c 

n 

E/k 

5a 

5.14E+  13 

0 

3464.02 

30a 

2.47E  +  28 

-3.37 

37644.69 

31a 

5.08E  +  23 

-2.51 

-34.02 

55a 

8.73E+  16 

-1.27 

-244.59 

mechanism;  thus,  a  single  reaction  involving  M  is  actually  several  reactions  with  the  same 
reaction  rate.  Table  3  gives  the  M  species  efficiencies  (a  value  of  unity  is  assumed  for  species 
not  listed).  These  factors  are  multiplied  by  the  concentration  of  each  species  in  the 
product-summation  terms  of  equation  6. 


Table  3.  Reaction  rate  efficiencies. 


Rxn 

M 

Efficiency 

Rxn 

M 

Efficiency 

19 

H20 

18.6 

31 

H20 

8.3 

19 

C02 

4.2 

31 

C02 

1.5 

19 

h2 

2.9 

41 

h2 

2.2 

19 

CO 

2.1 

41 

h2o 

6.7 

19 

n2 

1.3 

41 

co2 

3 

22 

h2o 

5 

41 

h2o 

5 

23 

h2 

0 

50 

n2 

1.48 

23 

co2 

0 

55 

02 

1.61 

23 

co2 

0 

55 

CO 

1.76 

30 

h2o 

4.4 

55 

co2 

4.17 

30 

co2 

2.3 

55 

h2o 

15 
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5.  Results  and  Discussion 


The  computational  domain  chosen  to  simulate  the  Penn  State  University  experiment  (figure  3)  is 
shown  in  figure  6  and  is  similar  to  that  used  to  simulate  the  ARL  experiment  sketched  in 
figure  2.  In  each  case,  this  domain  extends  from  the  capillary/extender-tube  on  the  left  to  the 
plate  on  the  right  (19  mm)  and  from  the  centerline  of  the  capillary  (and  plate)  to  a  fixed  radial 
distance  (40  mm)  that  is  determined  by  the  distance  from  the  plate  centerline  to  the  pressure 
probe  designated  PI  and  includes  a  small  radial  distance  beyond.  This  region  was  discretized 
using  1 54  axial  grid  cells  and  295  radial  grid  cells,  distributed  with  essentially  even  spacing 
throughout  as  shown  in  figure  6.  Some  degree  of  radial  grid  clustering  was  used  up  to  0.015  m 
in  order  to  more  accurately  resolve  the  formation  of  important  gas  dynamic  phenomena 
(expansions,  shocks  and  turbulent  mixing).  The  boundary  conditions  for  the  region  are 
symmetry  on  the  axis  (Y  =  0),  outflow  at  Y  =  0.04  m,  no-slip/no-penetration  on  the  capillary 
housing  surface  (X  =  0,  0.0015  <  Y  <  0.04  m)  and  on  the  plate  surface,  specified  inflow  at  the 
capillary  exit  (X  =  0,  0  <  Y  <0.015  m).  Initially,  the  entire  flowfield  is  filled  with  air  (0.8  mole 
fraction  of  N2  and  0.2  mole  fraction  of  O2).  Figure  7  shows  the  current  amplitude  variation  for  a 
typical  experiment.  Given  the  amplitude/time  and  the  physical  characteristics  of  the  capillary, 
the  capillary  model  (11)  generates  the  range  of  density,  velocity,  pressure,  temperature,  and 
species  distributions  shown  in  figures  8-10.  The  inflow  conditions  for  the  computational  domain 
follow  from  those  displayed  in  these  figures.  Figure  10  shows  that  there  are  five  prominent 
species  in  the  plasma  efflux  from  the  29  species  that  were  considered  in  the  capillary  model  (11). 
The  CFD  simulation  includes  39  species  (10  species  added  for  reaction  with  air):  electrons,  C, 

C+,  C++,  C“,  CH,  CH+,  CN,  CN+,  CO,  CO+,  C2,  C2+,  H,  H+,  I T’  H2,  H2+,  N,  N+,  NH,  NH+,  NO, 
NO+,  N2,  N2+,  O,  0+,  OH,  OH+,  02,  02+,  H20,  H02,  H2O2,  HNO2,  NO2,  C02,  and  03. 

To  aid  in  the  interpretation  of  computational  results  for  the  plasma  jet,  figure  1 1  shows  a 
schematic  of  the  gas  dynamic  features  expected  in  highly  underexpanded  supersonic  jets 
(27,  28).  The  efflux  of  plasma  from  the  inlet  generates  a  weak  precursor  shock  (A)  that  expands 
spherically.  Behind  this  shock  is  air;  the  plasma  is  entirely  contained  by  this  shock  and  is 
separated  from  the  air  by  an  irregularly  shaped  contact  surface  (B)  across  which  pressure  and 
velocity  are  preserved,  but  entropy  changes  discontinuously.  Expansion  waves  (Mach  cone) 
generate  at  the  inlet  (C),  travel  to  the  precursor  shock  (A),  are  reflected  as  weak  compression 
waves,  and  coalesce  into  a  strong  oblique  shock,  or  barrel  shock,  (D),  within  the  plasma  jet.  This 
barrel  shock  (D)  terminates  in  an  irregular  reflection  that  forms  a  triple-point  (E)  joining  the 
barrel  shock  (D),  its  reflection  (F),  and  a  normal  shock  (G)  or  Mach  disk.  Whereas,  the 
precursor  shock  (A)  is  relatively  weak  and  diffuse,  producing  a  mildly  supersonic  flow,  the 
barrel  shock  (D)  and  Mach  disk  (G)  are  strong  shocks  that  enclose  a  fully  supersonic  flow 
region. 
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Figure  6.  Computational  grid  used  for  the  simulation  of  the  Penn  State 

University  open-air  plate  experiment.  Lower  boundary  corresponds  to 
the  capillary  and  plate  centerline.  Partial  grid  shown  for  clarity. 


Figure  7.  Current  history. 
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Figure  8.  Efflux  density  and  velocity  histories. 


Figure  9.  Efflux  pressure  and  temperature  histories. 


Figures  12-18  show  the  time  history  of  the  plasma  jet  from  shortly  after  entrance  into  open-air 
(0.0208  ms)  to  long  after  impact  upon  the  plate  and  near  the  final  stages  of  efflux  from  the 
capillary  (0.241  ms).  Note  that  the  frames  in  figures  12-18  are  not  plotted  to  scale.  Mach 
number  contours  are  displayed  in  gray-scale  from  0  (bright  white)  to  1.5  and  above  (dark  black). 
Temperature  contours  are  displayed  in  gray-scale  from  300  K  and  below  (bright  white)  to 
28,000  K  and  above  (dark  black). 
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Figure  10.  Efflux  species  mole  fraction  histories. 


The  low-density  and  high-pressure  plasma  efflux  enters  the  open-air  as  a  highly  underexpanded 
jet.  An  expansion  wave  forms  at  the  capillary  exit.  A  precursor  shock  is  also  formed  as  the  jet  is 
moving  at  -6800  m/s  at  this  time,  or  -Mach  1 .7,  based  on  the  local  sound  speed  in  the  mixture 
(note  that  the  plasma  gas  is  of  very  low  molecular  weight,  -4.3  g/mole).  By  -0.0458  ms 
(figure  13),  the  predominant  gas  dynamic  features  in  the  jet  have  been  established,  including  the 
barrel  shocks  (seen  better  in  temperature  contours),  the  intersection  of  the  barrel  shock  with  the 
Mach  disk  (X  =  0.01 1  m,  Y  =  0.006  m,  figure  13a),  and  a  weak  reflected  shock  (seen  more 
clearly  in  figure  14a).  Due  to  the  variable  viscosity  in  the  flowfield  (see  equation  [5]),  the 
precursor  shock  is  more  diffuse  than  the  normal  shock  of  the  Mach  disk  or  the  oblique  barrel 
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Figure  12.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.0208  ms. 


Figure  13.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.0458  ms. 

shock  (note  the  subsonic  flow  region  between  the  normal  shock  and  the  precursor  shock  [light 
gray  region  in  Mach  number  and  dark  region  in  temperature  of  figure  13];  within  this  region 
resides  the  contact  surface).  The  precursor  shock  reaches  the  plate  at  -0.0458  ms.  Stagnation  of 
supersonic  flow  on  the  plate  causes  a  shock  reflection  that  moves  back  toward  the  capillary  as 
the  precursor  shock  travels  along  the  plate  (figures  13  and  14).  This  causes  the  Mach  disk  to 
recede  slightly  and  then  settle  at  a  position  0.01  m  from  the  capillary  (figure  14a).  At  this  stage, 
a  layer  of  high-pressure  and  high- temperature  gas  is  held  on  the  plate  near  the  centerline 
(i.e.,  plate  stagnation  region).  The  plasma  efflux  is  near  peak  output  (recall  figures  7-10). 

Figures  15-18  show  the  computed  results  for  times  when  the  plasma  efflux  is  waning; 
by  0.241  ms,  the  plasma  jet  has  ceased  (figure  18).  At  0.1 192  ms,  the  Mach  disk  is  greatly 
reduced  in  size  as  the  barrel  shock  dips  toward  the  centerline  and  the  reflected  shock  (at  the 
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Figure  14.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.0708  ms. 


(a) 


(b) 


Figure  15.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.1 192  ms. 

intersection  of  the  Mach  disk  and  the  barrel  shock)  impacts  the  plate  near  the  pressure-tap  P3 
(figure  15a).  Multiple  shock  reflections  from  the  plate  are  seen  to  subsequently  occur  (figure 
17a),  resulting  in  a  shroud  or  layer  of  high-pressure,  high-temperature  gas  along  the  plate 
surface.  As  the  plasma  efflux  ceases  (figure  18),  the  plate  slowly  begins  to  return  to  ambient 
conditions.  Figure  19  displays  these  events  as  time -history  plots  of  pressure  at  each  of  the 
pressure-tap  locations  (P1-P4).  An  additional  tap  location,  P0,  was  added  on  the  centerline  of 
the  plate,  but  has  no  experimental  equivalent.  The  precursor  shock  is  seen  to  arrive  at  the  plate 
(figure  19,  tap  P0)  at  -0.04  ms  and  subsequently  traverses  each  of  the  tap  locations  in  sequence. 
In  each  case,  the  pressure  is  momentarily  raised  to  10’s  or  100’s  of  MPa  and  then  rapidly  relaxes 
to  much  smaller  values.  As  further  flow  dynamics  proceed,  as  previously  described,  a  the  tap 
locations  re  again  pressurized  and  relaxed.  Eventually,  pressure  levels  at  all  tap  locations 
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Figure  16.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.1494  ms. 


Figure  17.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.1714  ms. 

relax  to  the  ambient  value  (not  shown  in  the  figure).  This  highly  transient  flow  has  important 
consequences  for  tests  of  the  ignition  and  sustained  combustion  of  a  solid  propellant  sample  to 
be  held  on  the  plate. 

Figures  20  and  21  show  comparisons  between  the  computed  and  measured  pressures.  Figure  20 
displays  the  computed  and  measured  pressure  peaks  (maximums)  and  the  corresponding 
time-of-arrival  of  the  precursor  shock.  The  results  for  a  19-mm  standoff  of  the  plate  compare 
well  with  the  measured  values.  In  general,  the  computed  shock  velocity  is  too  high,  perhaps  due 
to  the  absence  of  capillary  residue  (i.e.,  particles  from  the  vaporization  of  the  polyethylene  liner) 
in  the  simulation.  The  translation  of  these  particles  would  remove  energy  from  the  flow  and 
reduce  the  mean  velocity.  This  would  explain  the  overprediction  of  peak  pressures  at  the  P4 
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Figure  18.  Computed  (a)  Mach  number  contours  and  (b)  temperature  contours  at  time  =  0.2410  ms. 


Figure  19.  Computed  pressures  at  plate  taps. 

tap  location,  where  most  of  the  particles  are  usually  found  to  accumulate  (14).  Figure  21 
compares  the  P4  and  the  PI  tap  locations.  Predicted  results  for  the  P4  tap  location  have  been 
shifted  in  time  by  0.022  ms  in  order  to  line  up  the  arrival  time  for  the  precursor  shock.  With  this 
adjustment,  the  timing  of  subsequent  pressure  peaks  are  well  predicted  by  the  code,  but  the 
pressure  levels  are  too  low.  It  is  interesting  to  note  that  the  pressure  decrease  at  a  time  of 
0.17  ms  (figure  21a)  is  predicted.  This  pressure  drop  corresponds  to  a  flow  expansion  near  the 
P4  tap  location  at  this  time  (see  figure  17a).  The  results  for  the  PI  tap  location  are  also 
encouraging  (no  adjustments  have  been  made)  because  the  pressure  peaks  occur  at  regular  time 
intervals  for  both  the  predicted  and  measured  data;  however,  the  predicted  pressures  are  too  high. 
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Figure  20.  Peak  pressure  data  and  computations. 


Figure  21.  Computed  and  measured  pressure  histories  at  (a)  P4  and  (b)  PI  tap  locations. 


The  results  of  figures  20  and  21  demonstrate  that  the  CFD  code  is  able  to  track  the  dynamics  of 
the  flowfield  even  if  the  computed  pressure  levels  are  often  too  high.  Again,  the  absence  of 
particles  in  the  simulation  as  well  as  the  uncertainties  in  thermodynamic  and  transport  relations 
for  an  ionized  gas  at  these  temperatures  may  explain  some  of  the  discrepancies. 

Figures  22-24  display  some  of  the  chemical  aspects  of  the  flowfield.  Figure  22  shows  the  time 
histories  of  the  major  species  at  the  plate  centerline  (i.e.,  referred  to  the  P0  tap  location)  and  at 
the  tap  location  farthest  from  the  plate  centerline  (i.e.,  PI).  Comparison  of  these  figures  with 
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Figure  22.  Computed  time  variation  of  mole  fractions  at  two  locations  on  the  plate. 


Figure  23.  Computed  mole  fraction  contours  for  (a)  02  and  (b)  C  at  time  =  0.0262  ms. 


figure  10  (the  species  distributions  at  the  capillary  exit)  shows  that  a  large  proportion  of  the  ion 
species  are  actually  deposited  onto  the  plate  which  are  expected  to  ultimately  affect  the  ignition 
and  combustion  characteristics  of  propellant  mounted  on  the  plate.  Figure  23  and  figure  24a  are 
gray  scale  mole  fraction  contours  of  diatomic  oxygen,  atomic  carbon,  and  carbon  monoxide. 

These  results  show  that  while  the  plasma  is  devoid  of  O2  (large  white  region  extending  to  0.01  m 
in  figure  23  a),  there  is  some  O2  in  a  mixing  layer  on  the  edge  of  the  plasma  efflux  (note  region 
between  0.01  and  0.014  m)  as  well  as  atomic  carbon  that  has  diffused  outward  from  the  plasma 
jet.  These  combine  chemically  to  form  CO  (figure  24a)  in  this  mixing  region,  albeit  at  very  low 
levels  (1.0  x  10  10).  After  impact  with  the  plate,  CO  is  continually  formed  in  the  stagnation 
region  of  the  plate  (figure  24b).  It  can  be  concluded  that  some  mixing  between  the  plasma 


21 


Figure  24.  Computed  mole  fraction  contours  for  CO  at  times  =  (a)  0.0262  and  (b)  0.0458  ms. 

species  and  the  air  does  occur  and  that  this  mixing  promotes  a  low  level  of  chemical  reaction.  Of 
course,  these  results  are  dependent  on  the  plasma/air  chemical  mechanism  employed  in  this  study 
and  therefore  warrant  further  study. 


6.  Conclusions 


A  time-accurate  CFD  code  has  been  applied  to  the  modeling  of  the  high-temperature  chemically 
reactive  plasma  efflux  from  an  ETC  igniter  fired  into  open  air  and  impinging  on  an  instrumented 
plate.  The  CFD  code  has  been  linked  to  a  plasma  capillary  model  for  purposes  of  specification 
of  the  exit  conditions  of  the  igniter.  The  major  features  of  this  efflux  have  been  resolved  by 
numerical  simulation  revealing  a  highly  underexpanded  jet  with  a  strong  precursor  shock,  a 

barrel  shock  that  reflects  at  a  triple-point,  and  a  normal  shock  or  Mach  disk.  Impact  of  the  jet 
upon  a  plate  generates  a  stagnation  region,  a  reflected  shock  that  travels  back  toward  the 
capillary,  and  a  normal  shock  that  traverses  the  plate.  As  the  plasma  efflux  wanes,  the  barrel 
shock  reflects  from  the  centerline  and  then  the  plate,  causing  a  subsequent  series  of  pressure 
waves  along  the  plate.  Evidence  of  this  shock  system  and  pressure  wave  system  is  seen  in  both 
the  experimental  data  and  the  prediction.  The  presence  of  this  shock  system  has  important 
implications  for  mechanical  damage  to  solid  propellant  exposed  to  the  plasma  efflux.  Chemical 
conditions  at  the  plate  are  quite  different  from  those  at  the  ETC  igniter.  The  model  indicates  that 
ions  are  present  at  the  plate  and  that  mixing  (and  reaction)  between  the  plasma  species  and  the 
air  is  present  in  a  narrow  region  near  the  efflux  jet.  These  computer  simulations  represent  a 
unique,  detailed  modeling  of  ETC  igniter  plasma  efflux. 
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ADELPHI  MD  20783-1 197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  CS  IS  T 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 


1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
Only)  DTIC  OCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BEL  VOIR  VA  22060-6218 
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NO.  OF 

COPIES  ORGANIZATION 


ABERDEEN  PROVING  GROUND 


31  DIRUSARL 

AMSRD  ARL  WM  BD 
W  ANDERSON 
R  BEYER 
A  BRANT 
SBUNTE 
L  CHANG 
T  COFFEE 
J  COLBURN 
P  CONROY 
G  COOPER 
B  FORCH 
A  HORST 
S  HOWARD 
S  KARNA 
P  KASTE 
A  KOTLAR 
M  MCQUAID 
A  MIZIOLEK 
J  NEWILL 
M  NUSCA  (6  CPS) 

R  PESCE-RODRIQUEZ 
P  PLOSTINS 
BRICE 
JSAHU 
S  SILTON 
P  WEINACHT 
A  WILLIAMS 


27 


Intentionally  left  blank. 


28 


